arXiv:1504.07630vl [cond-mat.mes-hall] 28 Apr 2015 


Spatially resolved edge currents and guided-wave electronic 

states in graphene 


M. T. Allen\ O. Shtanko^, 1. C. Fulga^, A. Akhmerov^, K. Watanabi^, T. Taniguchi^, P. 
Jarillo-Herrero^, L. S. Levitov^, and A. Yacoby^* 


^ Department of Physics, Harvard University. Cambridge, MA 02138. 

^ Department of Physics, Massachusetts Institute of Technology. Cambridge, MA 02139. 

^ Department of Condensed Matter Physics, Weizmann Institute of Science. Rehovot, 

Israel. 

^ Kavli Institute of Nanoscience. Delft University of Technology. Lorentzweg 1. 2628 CJ 

Delft. The Netherlands. 

^ Environment and Energy Materials Division, National Institute for Materials Science. 1-1 

Namiki, Tsukuba, Ibaraki, 305-0044 Japan. 


* Corresponding author. E-mail: yacoby@physics.harvard.edu 



1 


Abstract: A far-reaching goal of graphene research is exploiting the unique properties of 
carriers to realize extreme nonclassical electronic transport. Of particular interest is harness¬ 
ing wavelike carriers to guide and direct them on submicron scales, similar to light in optical 
fibers. Such modes, while long anticipated, have never been demonstrated experimentally. 

In order to explore this behavior, we employ superconducting interferometry in a graphene 
Josephs on junction to reconstruct the real-space supercurrent density using Fourier methods. 

Our measurements reveal charge flow guided along crystal boundaries close to charge neu¬ 
trality. We interpret the observed edge currents in terms of guided-wave states, confined to 
the edge by band bending and transmitted as plane waves. As a direct analog of refraction- 
based confinement of light in optical fibers, such nonclassical states afford new means for 
information transduction and processing at the nanoscale. 

Electrons in Dirac materials such as graphene can be manipulated using external fields to control 
electron refraction and transmission in the same way that optical interfaces in mirrors and lenses can 
manipulate light [1^]. Several of the key ingredients, including phase-coherent Klein transmission 
and reflection [5-7], ballistic transport [8] and transverse focusing on micrometer scales [9], have 
already been established. One promising yet unexplored direction, which we investigate here, is the 
quasi-ID confinement of electrons in direct analogy to refraction-based confinement of photons in 
optical fibers. Electronic guided modes formed by a line gate potential, while discussed in the litera¬ 
ture [10-13], have so far evaded direct experimental detection. Extending the fiber optics techniques 
to the electronic domain is key for achieving control of electron waves at a level comparable to that 
for light in optical communication systems. 

Rather than pursuing the schemes discussed in Refs. [10-13], here we explore modes at the 
graphene edges. The atomically sharp graphene edges provide a natural vehicle for band bending 
near the boundary which then confines the electronic waves in the direction transverse to the edge. 
The resulting guided “fiber-optic” modes propagate along the edge as plane waves, decaying into the 
graphene bulk as evanescent waves. In analogy with optical fibers, such modes are situated outside 
the Dirac continuum (see Eig. 1 A,B). As discussed below, the mode frequency in monolayer graphene 
is 

(1) uj = v\k\ — |i;| <-u = 10® m/s, 

where the damping 7(A:) accounts for scattering by disorder. Eor the practically interesting regime of 
disorder originating from edge roughness, the damping is expected to quickly vanish at long electron 
wavelengths near charge neutrality, scaling as 7(A:) ~ A^. In this regime, as discussed below, mode 
scattering by the edge is near-specular and is thus immune to backscattering. This approach to car¬ 
rier guiding is particularly appealing because of the ease with which band bending at the graphene 
edge can be realized, as well as because there is no threshold for fiber-optic states to occur: they are 
induced by an edge potential of either sign, positive or negative, no matter how weak (see discussion 
below and in the Supplementary Methods). The presence of such guided modes enhances the density 
of current-carrying states at the edge for doping near charge neutrality, while uniform behavior is 
recovered at higher carrier concentrations (see Eig. 1C and Eig. SI). 

The edge currents associated with guided states, anticipated at zero magnetic field, have so far 
eluded experimental detection due to the challenge of imaging current with submicron spatial res¬ 
olution. In particular, scanning tunneling spectroscopy (STS) images density of states but not cur¬ 
rent flow [14,15], while macroscopic conductivity cannot distinguish the edge and bulk contribu¬ 
tions [16, 17]. Additionally, these techniques covered disparate length scales, with STS probing 
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atomic-scale lengths and transport covering micron to millimeter distances. High resolution den¬ 
sity of states measurements along unzipped carbon nanotube boundaries and graphene quantum dots 
were obtained spectroscopically using STS, but the dispersive nature and current-carrying capacity of 
these states was unresolved [14,15]. 

With this motivation, we developed a technique to spatially image electric current pathways and 
applied it to high-mobility graphene. We employ Fraunhofer interferometry in a graphene Josephson 
junction to reconstruct the spatial structure of the electronic states which transmit supercurrent. To 
implement this approach, we measure gated Josephson junctions consisting of graphene coupled to 
superconducting titanium/aluminum electrodes (Fig. ID). A gate electrode is used to tune the carrier 
density n. In order to access the intrinsic properties of graphene at densities near charge neutrality, we 
isolate the flakes from substrate-induced disorder through placement on thin hexagonal boron nitride 
(hBN) substrates [18]. We study four bilayer and one monolayer devices, all of which exhibit similar 
behavior (Table SI). Figures lE-H exemplify transport data from one of the bilayer devices described 
above. Upon sweeping DC current bias Idc^ a sharp transition in resistance between dissipationless 
and normal metal behavior appears at a critical current Ic, a signature of the Josephson effect (Fig. 
1E,F). 

We employ superconducting quantum interference to extract a spatially resolved image of the su¬ 
percurrent density across the flake. It is possible to obtain real-space information because application 
of a magnetic flux <h through the junction area induces a position-dependent superconducting phase 
difference A(j){x) = 2'k^x/^qW parallel to the graphene/contact interface [19], where 4>o = /i/2e 
is the flux quantum, h is Planck’s constant, e is the elementary charge, and W is the width of the 
flake (Eig. ID). The resulting interference is displayed in Eigure IE, a plot of differential resistance 
dV/dl as a function of Idc and magnetic field B. Measurements of the AC voltage drop dV across 
the junction in response to an AC current modulation dl were conducted using lockin techniques 
in a dilution refrigerator at 10 mK, well below the critical temperature of Al. The critical current Ic, 
obtained by extracting the value of loc at the maximum derivative dV/dl, exhibits modulations in B 
field that arise from Eraunhofer-like magnetic interference. As the flux threading the junction winds 
the superconducting phase along the length of the contact, the critical current Ic can be expressed 
quantitatively as the magnitude of the complex Eourier transform of the current density distribution 
J(x). That is, Ic = \Ic{B)\, where 

I-WI2 

(2) Ic{B) = / J(x) • 

J-WI2 

where L is the distance between contacts. Relevant for wide junctions {L <C W) such that the cur¬ 
rent is only a function of one coordinate. Equation (2) provides a simple and concise description of 
our system. The spatial distribution of supercurrent thus dictates the shape of the interference pat¬ 
tern [19-21]. 

Our results, obtained with this technique, show a strikingly different behavior at high and low car¬ 
rier densities. We observe conventional Josephson behavior with uniform current flow at high density, 
for which the normalized critical current Ic{B)/lc{0) ~ | sin(7r<I>/<ho)/(vr<h/<ho)| is described by 
single-slit Eraunhofer diffraction (Eig. IE). Defining features of such interference include a central 
lobe of width 24>o and side lobes with period <ho and decaying 1/B amplitude. However, near the 
Dirac point, our results exhibit a striking departure from this picture and show an enhancd “SQUID- 
like” interference (Eig. IE) [22]. Such behavior arises when supercurrent is confined to edge channels 
and is characterized by slowly decaying sinusoidal oscillations of period 4>o. Importantly, these two 
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regimes are easily distinguishable without mueh analysis by observing the width of the eentral lobe 
whieh is twiee as wide for the uniform ease as eompared to the ease of edge flow. 

The real-spaee eurrent distribution ean be quantitatively obtained by inverting the relation in Eq.(l) 
with the help of the Fourier teehniques of Dynes and Fulton [20] (see Supplementary Methods). A 
more numerieally expensive Bayesian estimation of the eurrent distribution produees eurrent distribu¬ 
tions and standard error estimates that agree with the Fourier teehniques (see Supplementary Methods 
and Fig. S2). The resulting eleetron density map reveals strong eonfinement of supereurrent to the 
edges of the erystal near the Dirae point (Fig. IH). This phenomena is a robust experimental feature 
seen in all five deviees. The width of the edge eurrents, extraeted quantitatively from Gaussian fits, is 
on the order of the eleetron wavelength (~ 200 nm) and roughly eonsistent aeross multiple samples. 
This value is likely an upper bound beeause the peak width is manifested in the deeay envelope of the 
interferenee pattern, external faetors that suppress the eritieal eurrent amplitude at high B, sueh as ae- 
tivation and deeay of the A1 supereondueting gap, may also eontribute to peak broadening. Although 
the maximum fields of B ~ 10 — 25 mT are a small fraetion of the eritieal field of A1 {He ~ 100 mT), 
the aetual edge states may be narrower. At high eleetron density, eonventional single-slit Fraunhofer 
interferenee is reeovered (Fig. IE), suggesting a uniform distribution of supereurrent (Fig. IG). 

By tuning earrier eoneentration with a gate eleetrode, our measurements reveal eoexistanee of edge 
and bulk modes at intermediate densities in (Fig. 2, monolayer graphene). Starting at the Dirae point, 
the deviee exhibits SQUID-like quantum interferenee through eharge neutrality and the spatial image 
of supereurrent reveals edge-dominated transport (Fig. 2A-C). As density is inereased, bulk eurrent 
flow inereases, erossing over to mostly uniform flow (Fig. 2D) and eonventional Fraunhofer interfer¬ 
enee at high eleetron density (Fig. 2E). To traek the evolution of edge and bulk eurrents with density, 
we plot line euts of the individual eontributions in Figure 2E, where the gate voltage eorresponding 
to the charge neutrality point is identified as a dip in the current amplitude. Notably, raw interference 
near the Dirac point (Fig. 2E) and at high electron concentration (Fig. 2F) exhibit the salient features 
that distinguish edge-dominated from bulk-dominated transport, including a width of (j)Q versus 2(j)Q 
of the central lobe, as well as Gaussian versus 1/B decay of the lobe amplitudes for low and high 
densities, respectively. Finally, due to the spin and momentum conserving nature of Andreev reflec¬ 
tion, the observed supercurrents suggest an absence of strong spin and valley scattering at the edge in 
high quality samples. 

Similarly, we systematically explore the correspondence between edge and bulk flow in bilayer 
graphene and detect boundary currents in the presence of broken crystal inversion symmetry (Fig. 
3). As the Fermi energy approaches the Dirac point from the hole side, the bulk is suppressed faster 
than the edge, leading to emergence of robust edge currents near zero carrier density. Upon tuning 
the Fermi level to high electron doping, a uniform distribution re-emerges (Fig. 3A,B). We note that 
the range in hole density over which the bulk contribution is recovered varies in different devices. 
At hole doping, an intrinsic p-n junction emerges at the graphene-electrode boundary due to contact- 
induced charge transfer, a phenomenon intensively studied using both transport and photocurrent 
techniques [23-25]. In this regime, we find that edge transmission is augmented relative to the bulk, 
perhaps because an intrinsic gap at the p-n interface can block bulk flow more effectively in bilayer 
graphene. (In this device, current distributions are not plotted at the immediate Dirac point due to sup¬ 
pression of proximity-induced superconductivity at high normal state resistances.) Furthermore, we 
apply an interlayer electric field E to break crystal inversion symmetry and induce a bandgap [26,27], 
manifested as a gate-tunable insulating state at the Dirac point (Fig. 3C,D). In this regime, we find 
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that conductance is mediated by edge currents that enclose the bulk, even in the presence of a field- 
induced gap (Fig. 3E,F). 

Our measurements establish that edge currents dominate electronic transport in graphene at zero 
magnetic field near the Dirac point. As demonstrated below, such edge currents can arise from ‘fiber¬ 
optic’ electronic states formed in a line potential localized near the edge. As a simple model exhibiting 
an analogy between electronic guided states and refraction-based confinement of light in fiber optics, 
we consider massless Dirac particles in graphene monolayer in the presence of a line potential: 

(3) e^{x,y) = {vap+ V{x)) uRilO^m/s. 

We seek plane-wave solutions of the form Tp{x,y) = Here k is the wavevector compo¬ 

nent along the line and (j){x) is a two-component spinor wavefunction depending on the transverse 
coordinate, and ai are Pauli matrices. This problem can be tackled by a gauge transformation 

(4) (^(x) = 6(x) = ^ [ V{x')dx'. 

it-V Jq 

Such a transformation generates a mass term in the Dirac equation and eliminates the potential V (x) 
entirely. We further simplify the resulting equation using the identity = ay cos(20) — 

az sin(20) to obtain e^{x) = hv{—iaxdx + kay cos(26*(y)) — kaz sm{29{x))(f>{x). 

As a simple example, we consider the case of an armchair edge, for which the problem on a half 
plane for carriers in valleys K and K' is equivalent to the problem on a full plane for a single valley. 
Applying the above equation to a potential localized in an interval —d < x < d and focusing on the 
long-wavelength modes such that kd <C 1, we can approximate d{x) ^ sgn (x) with the parameter 

u = -^ f^^V(x')dx'. We arrive at the seminal Jackiw-Rebbi problem for a Dirac equation with a 
mass kink 

(5) e^(x) = hv(—iaxdx-h ayk + azm(x))^(x), k = kcosu, m(x) = —ksmusgn(x). 

The Jackiw-Rebbi problem can be solved exactly and explicitly [28]. Guided-wave states for such a 
Hamiltonian are described as products of the zero-mode state found for fc = 0 and the plane wave 
factors The energies of these states are simply e = ^hvk with the sign given by sgn (m(0-|-) — 
m(0—)). This gives the dispersion in Eq.(l) with v = ±vcosu and the sign sgn (sin tt). Since 
|x| < X, for each k the energies of these states lie outside the bulk continuum |e| > hv\k\, which 
ensures decoupling from the bulk states and confinement to the region near the x = 0 line. The 
connection with the theory of zero modes indicates the robustness of such confinement. Similar fiber¬ 
optic states are obtained for an edge potential in graphene bilayer (see Eig. IB). 

To assess compatibility with this model, we compare supercurrent density measurements with a 
theoretical prediction of density of states (Eig. 4). Real space line cuts of current flow J(x) in bilayer 
device BL3 at fixed densities are provided in Eig. 4a, showing edge currents near the Dirac point and 
a continuous evolution of bulk flow. Theoretical plots of density of states as a function of position 
(Eig. 4B), obtained from the ‘fiber-optics’ model of transport, exhibit qualitatively similar behavior. 
Eor the simulation, an effective delta function potential approximation is used with the best-fit value 
A = 0.5 eV-nm (see Supplementary Methods). 

Despite the edge roughness inevitable in our devices, the observed edge currents are robust. This is 
consistent with the guided-wave model. Indeed, the long-wavelength guided modes are characterized 
by a large transverse confinement lengthscale, since the evanescent wave decay defines a transverse 
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confinement lengthscale on the order of electron wavelength A, which is much larger than the edge 
roughness scale In this case, most of the mode wavefunction resides far from the edge, at dis¬ 
tances r ~ A and is thus not susceptible to disorder scattering. For typical wavelength values 

A ~ 10 — 100 nm, which exceeds the atomic scale characteristic for the edge disorder by two to three 
orders of magnitude, we expect backscattering suppression by as much as (A/^)^ ~ 10^ — 10® times 
to be readily achievable [see Eq.(l) and Supplementary Methods]. This is similar to optical guided 
states in so-called “low-contrast” optical fibers, which feature a similar suppression of disorder scat¬ 
tering by a large ‘Thomson factor’ A/^. This suppression is key for achieving extremely long mean 
free paths in such fibers. 

To the best of our knowledge, the fiber-optic modes is the only model which is fully consistent 
with the observations. For instance, edge density accumulation can influence the edge potential and 
in principle also support guided edge currents. However, the fact that the charge neutrality points for 
both edge and bulk roughly coincide in density n suggests an absence of a positional charge imbal¬ 
ance on a large scale (Fig. 2F). In addition, edge-dominated current flow is observed near the Dirac 
point but not at higher densities, the behavior not expected for strong edge doping. Arguments against 
random density inhomogeneities (of the kind familiar from our previous studies of electron-hole pud¬ 
dles) include the reproducibility of edge currents with width on the order of electron wavelength 
across many samples, as well as the observation that edge currents tend to be stronger in clean sam¬ 
ples with ballistic Fabry-Perot interference. Large charge inhomogeneities across the sample would 
suppress Fabry-Perot interference and are thus unlikely. 

One appealing aspect of the fiber-optic model is that it can naturally accommodate a wide range of 
different microscopic physical mechanisms discussed theoretically in the literature [28-32] that may 
produce an edge potential. Examples include pinning of the Eermi energy to the low-energy states 
due to broken A/B sublattice symmetry [30-32], density accumulation caused by dangling bonds or 
trapped charges at the boundaries, or electrostatics [28,29]. The competition of these effects can 
produce a complex dependence of the edge potential V(x) on carrier density. Pinpointing the precise 
microscopic origins of the edge potential requires further study. 

Lastly, it is widely known that the A/B sublattice imbalance for broken bonds at the edge can lead 
to edge modes in pristine graphene at neutrality. Such “dispersing zero-mode states” can exist even 
in the absence of a line potential, forming edge modes for an atomically perfect zigzag edge [30-32]. 
However, our simulations for disordered edge show that these states are highly localized on the dis¬ 
order length scale, and also that edge roughness induces strong scattering between the states at the 
boundary and in the bulk, producing relatively short mean free path values and hindering ballistic 
propagation. We therefore conclude that such states are unlikely to contribute to the observed bound¬ 
ary currents. 

To summarize, we present evidence of zero field edge currents in a graphene Josephson junction. 
These findings underscore the relevance of edge effects for graphene transport and provide a new tool 
for exploration of boundary currents, a topic of emerging interest. Eormed due to band bending at the 
edge, the observed guided states demonstrate confinement of electron waves at a level comparable 
to that for light in photonic systems. This defines a new mode for transmission of electronic signals 
at the nanoscale. We also anticipate this work will inspire more detailed investigations of boundary 
states in graphene and other materials. Our techniques open the door to fast spatial imaging of current 
distributions along more complicated networks of domains in larger crystals, increasingly relevant as 
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graphene sheets are sealed to larger dimensions [33]. Sueh eapahilities are also of fundamental inter¬ 
est due to the predieted topologieal nature of edge states along staeking domain boundaries in bilayer 
graphene [34,35]. 
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Figure legends 

Fig. 1. ‘Fiber-optic’ modes and spatially resolved current imaging in a graphene Josephson 
junction. (A, B) Guided edge modes indueed by an intrinsie band bending near erystal boundary, 
for single-layer and bilayer graphene (sehematie). Mode frequeneies positioned outside the Dirae 
eontinuum ensure mode deeoupling from the bulk states. Guided modes exist for any edge potential 
no matter how weak. In a single layer, mode veloeity ehanges sign as the potential strength inereases, 
see Eq.(5). In a bilayer, the modes oeeur in pairs [green and red curves: dispersion for positive and 
negative potential strength, respeetively]. (C) The guided modes are manifested through peaks in 
the density of eurrent-earrying states at the erystal boundaries, prominent near eharge neutrality (red: 
n = 0.05 X lO^^crn”^; blue: n = 2.5 x lO^^crn”^). (D) Sehematies of supereondueting interferom¬ 
etry in a graphene Josephson junetion, whieh is used to image the spatial strueture of eurrent-earrying 
states. A flux is threaded through the junetion area to produee interferenee patterns, as eurrent bias Vsd 
is applied through the supereondueting eleetrodes and the voltage drop aeross the deviee is reeorded. 
Carrier density n is tuned by a gate voltage kb- (E, F) The reeorded interferenee pattern is of a 
single-slit Fraunhofer type at high earrier density, turning into a SQUID-like interferenee near neu¬ 
trality (eolorseale is dV/dl (fJ) for deviee BLl). (G, H) Current flow, extraeted from the interferenee 
data using Fourier teehniques, is uniform at high earrier density and peaks at the erystal edges for ear¬ 
rier density elose to neutrality. 

Fig. 2. Gate-tunable evolution of edge and bulk current-carrying states in graphene. (A) 

Edge-dominated SQUID-like interferenee pattern at neutrality in deviee MLl (n = 2.38 x 10® em~^; 
eolorseale is dV/dI{Q)). (B, C) Real-spaee image of eurrent flow eonfined to the boundaries over a 
range of densities near neutrality, shown alongside with the raw interferenee data (eorresponding to 
the white box in (D)). (D) A real-spaee map of eurrent flow as a funetion of eleetron eoneentration 
reveals eoexistenee of edge and bulk modes at intermediate densities. (E) Conventional Eraunhofer 
pattern for uniform eurrent flow at high eleetron density (n = 7 x 10^^ em“^). (F) Comparison of 
eurrent amplitudes along the edge (red) and bulk (blue) from the plot in panel (C). Current flow is 
edge-dominated near neutrality. Note that minima for both eontributions eoineide in n, indieating 
that a positional edge/bulk density offset is not present. 

Fig. 3. Boundary currents in bilayer graphene in the presence of broken crystal inver¬ 
sion symmetry. (A) Spatially resolved supereurrent map in deviee BL2, in a normalized plot of 
J{x)/Jmax{x)- Edge-dominated transport oeeurs near eharge neutrality, while an inereasing bulk 
eontribution is tuned with earrier eoneentration. (B) Comparison of eurrent amplitudes along the edge 
(red) and through the bulk (blue) from panel (A). Enhaneed edge eurrents are prominent at neutrality, 
whereas a uniformly distributed flow is reeovered at high densities. The normal state eonduetanee 
G{e^/h) vs. earrier density is also shown (blaek). (C) Measurement sehematie for supereondueting 
interferometry in a dual-gated bilayer graphene Josephson junetion. A dual-gated deviee eonsists of 
bilayer graphene flake on hBN with a suspended top gate, where applieation of voltages Vt and Ub 
on the top and baek gates enables independent eontrol of the transverse eleetrie field E and earrier 
density n. (D) Resistanee map as a funetion of Ft, and Vt for bilayer BL4. Enhaneed resistanee at 
high E fields indieates the emergenee of a gate-tunable insulating state due to broken erystal inversion 
symmetry. (E) Spatially-resolved boundary eurrents as a funetion of E field. The vertieal axis is a 
traee along the red path labeled in (B). (F) Sequenee of Fraunhofer measurements at various loeations 
on the eurrent map in panel (E). 



Fig. 4. ‘Fiber-optics’ theoretical model of transport in graphene. (A) Real-space maps of 
measured current flow J{x) in bilayer device BL3 at fixed carrier densities on the hole side, showing 
edge currents near the Dirac point and a continuous evolution towards bulk flow. (B) Theoretical plot 
of spatially resolved density of states in bilayer graphene at fixed carrier densities for edge waveg¬ 
uide model. For the simulation, an effective delta function potential approximation is used with the 
best-fit value A = 0.5 eV-nm (see SOM). Band mass of bilayer graphene is taken 0.04me where rUe 
is electron mass. 
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Figure 4 
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Materials and Methods 

Modeling electronic guided modes 


A full model of supercurrent-carrying states in our system should account, in principle, 
for a number of microscopic effects. This includes, in particular, the microscopic details of 
transport through the NS interfaces, the realistic edge potential profile due to band bending 
near graphene edge, as well as the effects of disorder. Since treating all these issues simulta¬ 
neously and on equal footing makes such a modeling a daunting task, here we resort to some 
simplifications. First, we will completely ignore the effects of induced superconductivity, fo¬ 
cusing on the normal metallic state of a pure graphene. Second, we consider a clean system 
and account for disorder scattering perturbatively at the end. Third, since states in a clean 
system, being delocalized, are capable of carrying supercurrent, we will focus on evaluating 
the density of states (DOS) taking it to reflect on the current-carrying capacity of the system. 
Of course, such an approach may be questioned for disordered systems in which some states 
are localized, and therefore can contribute to DOS but not to supercurrent. However, taking 
into account that in a clean system all states possess a roughly similar current-carrying ca¬ 
pacity, we adopt this approximation on the merit of its simplicity. 

Turning to the discussion of system geometry, we note two points. First, as discussed in the 
main text, the problem of guided states on a halfplane near the edge x > xq can be mapped 
onto a similar problem on a full plane by accounting for the states in valleys K and K' mix¬ 
ing at the edge. This mapping is particularly transparent for the armchair edge, where the 
boundary condition for the spinor wavefunctions in the two valleys is simply ifjx + V’ir' = 0. 
In this case, one can see that the two-valley half-plane problem is mathematically equiva¬ 
lent to the problem posed on a full plane for particles in just one valley, provided the line 
potential for the latter problem is taken to be a sum of the original edge potential and its 
mirror-reflected double, V{x > Xq) —)■ V{\x — XqI). 

Second, the states with the wavelengths larger than the edge potential width can be de¬ 
scribed by a delta function approximation. In that, a realistic microscopic potential V(x) is 
replaced by a delta-function pseudopotential V{x) = \6{x — xq), where A = / V{x')dx' 
and Xq is the edge position. For a system of width w with two parallel edges positioned at 
Xq = ±w/2 we therefore arrive at the model 

(1) V{x) = X5{x + w/2) + X5{x — w/2), 

with —oo < X < oo. Carriers in this system are described by the massless Dirac Hamiltonian 

(2) H = Hq + V (x), Hq = vdip^ + va2Py 

with V ~ 10®m/s the carrier velocity and ai ^2 the pseudospin Pauli matrices. 

As stated above, we will use spatially-resolved DOS for the problem (2) as a measure 
of current-carrying capacity of the system. In justification we note that an electron system 
carrying normal electric current can be understood in terms of changes in the occupancy of 
the states near the Fermi energy. As a result, the spatially-resolved current density will vary 
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in the same manner as DOS 

(3) iV(/i,r) = n(r) = (V^^'(r)^(r)). 

dfi 

Here n is the total earner density and /i is chemical potential. Below we evaluate DOS as 
a function of position and energy, focusing on the characteristic features due to the guided 
modes. 


Taking into account that typical wavelength values of relevant electronic states, A ~ lOnm 
, are much smaller than the distance between edges w ~ l/im, we can represent DOS in the 
form 


(4) = iVo(/r) + -w/2) + Ni{fi,x + w/2) 


where Nq is the DOS of a uniform infinite system. 


(5) 


N.(e) = 


2'Kh?v‘^ 

and Ni is the contribution to DOS from a single delta-function line potential, placed at a; = 0. 
Below we derive an expression 


(6) Ni{e,x) = 


4A 


Im 


dp 




■nhv J 2tI Ke,p [4A£ -f (4 - X‘^)hvKe^p] 


, Ke,p = y/p2 - {e/hvy 


where the energy e is taken to have an infinitesimal positive imaginary part. In the final 
result for DOS e must be replaced by the chemical potential, e = p. The spatial dependence 
described by Eq.(6) is shown on Fig. SI. 


In our model, which is essentially non-interacting, the effects of screening can be included 
ad hoc by treating the potential strength in Eq.(l) as a function of carrier density. Since the 
latter is parameterized by /i, we will use a simple model 

C7) A —^ A = --^ 

1 -f (|/i|/po)" 

where the parameter po depends on microscopic details. Comparing to the data indicates 
that a reasonably good fit can be achieved for a ^ 2. 


Modeling results are presented in Fig. 1(c) of the main text for energies corresponding to 
carrier densities n = 0.05 ■ lO^^cm”^ (red curve) and n = 2.5 ■ lO^^cm”^ (blue curve), 
where we evaluated n accounting for the spin and valley degeneracy in graphene. Poten¬ 
tial strength is chosen to be A = —1.5 hv ^ 1 eV-nm and the screening parameter value is 
Pq = 0.2y'Tih?v'^no ~ 7 meV, where Uq = lO^^cm”^ is the corresponding scale for density. 


The simulation for graphene bilayer, which was used to generate Fig. 1(b) and Fig. 4(b) of 
the main text, was carried out using an effective delta function potential approximation, as 
above. Greens function expressed through the T-matrix was used to obtain mode dispersion 
and DOS in a manner similar to our treatment of modes in a single layer. For the delta 
function strength we used the best-fit value A = 0.5 eV-nm (and no screening). 
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Microscopic derivation 


Here we consider long-wavelength modes for a potential line positioned at x = 0. This 
problem is described by the Hamiltonian (2) with V(x) = XS{x). For this problem we con¬ 
struct the Greens function which takes the full account of scattering by the potential. As 
is well known, the discrete spectrum of the system (in our case, the guided modes) can be 
conveniently expressed through the poles of the electron Greens function. Likewise, the 
spatially-resolved DOS is expressed as the Greens function trace. The Greens function, in 
turn, can be straightforwardly evaluated using Dysons’s equation and the T-matrix represen¬ 
tation: 


G — Gq -f GqVGq -f GqVGqVGq -!-■■■ — Gq -I- GqTGq 


( 8 ) 


where Go = {is — Hq)~^. 

Assuming that the phase and amplitude of the electron wavefunction are given by a con¬ 
tinuous function of x, we can express the quantity T as 



(9) 


The continuity assumption should in practice be relaxed by a weaker assumption accounting 
for the phase jump of the wavefunction across the delta function potential at x = 0 (see main 
text). Here, however, we will proceed with Eq.(9) on the merit of its simplicity. Evaluating 
the integral in Eq.(9) gives 



( 10 ) 


where we defined 


hvPy 


( 11 ) 




-f 


Here e is the Matsubara frequency, with a suitable analytic continuation ie ^ e + it) to be 
performed at the end. 

The T-matrix poles give the guided modes dispersion 


-f 


e = Aihu\py\, u = V 


( 12 ) 


where the sign is given by ± = signA. Since |m| < v, the energies £ = ±u\py \ are positioned, 
for each py value, outside the Dirac continuum of the bulk states. This expression behaves 
in a qualitatively similar way to the exact dispersion derived in the main text, Eq.(l) [see 
Eig.l(a) of the main text]. The guided modes described by Eq.(12) are quasi-lD states that 
propagate as plane waves in the y direction along the x = 0 line and decay exponentially as 
evanescent waves in the transverse direction. 

Spatially-resolved DOS can be evaluated as 


N{e^ r) 


1 


ImTr G(£, r, r')r=r' 


TT 


(13) 
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where the energy variable is analytieally eontinued from positive imaginary to real values via 
ie —)■ £ + iO and a trace is taken over pseudospin variables. To proceed with our calculation, 
we will need Greens function evaluated in a mixed position-momentum representation 

(14) Go(£,Pj,i) = y ^e'’’‘”Go(£,p) 

where K{ie) = ^^{e/hvy + 

The trace of an equal-point Greens function in Eq.(13) then could be evaluated from Eq.(8) 
with the help of Eq.(lO): 


(15) 


i:iG{e,x' = x) = 

Pv 



AXp^e 


2 -2K\x\/h 


hv [(2 -I- iXey — 


where the two terms represent contributions of Go and GqVG o, respectively. 


As a warmup, we consider the first term of (15). Introducing a UV cutoff po = £o/^v we 
evaluate the sum over py as 

dpy e 


(16) 


In^. 


/_pg 27rh -f Tihv e 

Performing analytic continuation e ^ 5 — ie, we arrive at 

^0 


(17) 


Noe) = 


Tm In 


5 — is 

where 5 = -fO. Taking the imaginary part, we obtain the expression in Eq.(5). 


Next, we proceed to evaluate the second term in Eq.(15). Performing the same analytic 
continuation, we arrive at the result in Eq.(6). The expression in Eq.(6) can be conveniently 
analyzed by dividing the integral into two parts, taken over the domains \py\ > e/hv and 
\py\ < e/hv, respectively. The latter contribution is particularly simple because it is governed 
by the pole (12) and can be easily evaluated, giving 

ng'i N (ft)— 

^ h?vu{A-xY 


This contribution is solely due to the guided edge mode. As illustrated in the Pig SI, this 
term dominates the peak structure in DOS for guided waves. 


We used the full expression in Eq.(6) to produce the spatially-resolved DOS curves shown 
in Pig. 1(c) of the main text. In that, we accounted for screening, as described in Eq.(7). 
Because of screening, the peak structure is more prominent at low chemical potential, and is 
suppressed relatively to the bulk DOS at high chemical potential values. 


The effect of disorder 













6 


Here we estimate the disorder scattering rate 7 (A;) for guided modes [see Eq.(l) in the 
main text and accompanying discussion]. We will model edge roughness by a fluctuating 
delta function strength, treating the fluctuations as a gaussian white noise: 

(19) l/(a;,|/) = (A + (5A(|/)) (5(a;), {5\{y)5\{y')) = a5{y - y). 

Writing the Greens function as a series in the potential V + SV, Eq.(19), we have 

(20) G = Go + Go(V + 6V)Go + Gq(V + SV)Go(V + SV)Go + ... 

Averaging the Greens function over disorder, we only need to account for the pair correlators 
{6X{y)6\{y')). In a non-crossing approximation, we express the disorder-averaged Greens 
function through a suitable self-energy 

(21) {G) = Go + Go(y -\- S)Go -f GoiV + S)Go(E -f S)Go -f ... 

where 

( 22 ) J:{e) = a j ^G{e,py,x,x')^=^>=o 

The quantity (22) is complex-valued, with the imaginary part expressed through the density 
of states at a: = 0 as 

(23) ImS(e) = —7raN{e)x=o 

The disorder scattering rate for the guided waves can now be found from the dispersion 
relation obtained from the T-matrix pole, Eg(9), which is corrected by the presence of S as 
follows 

(24) 1 + (A + E(*£))?^J^ = 0. 

Here we continue to use Matsubara notation, as in Eqs.(9),(10). 

Since the density of states scales linearly with energy, N{e) ~ |e|, we can solve Eq.(24) in 
the long-wavelength limit treating S(ie) as a perturbation. Writing e = eo{py) + Se, where 
^0 = u\py \ is a solution for S = 0, we linearize in 5e to obtain 

1 / \ 

(25) 6e =-—jJ:{ieo)\py\ 

After analytic continuation, we obtain 

(26) 7(P?;) = j^(l- 

Accounting for the linear scaling N ~ |£|, we find that the damping rate scales as a square 

of Py, 

( 27 ) X{Py) = 

at small py. A similar dependence, albeit with a different prefactor, is found at large Py. 

Erom this we conclude that the modes are undamped over lengthscales ~ A^/.^, where A 
is a wavelength and is a disorder lengthscale. Taking realistic values A ~ 10 — 100 nm 
and ^ 0.1 nm, we obtain an estimate for the guided mode mean free path in the 1 — 10 /rm 
range. These large values can be traced to the weak confinement of the waves at small Py. 
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The weak confinement results in the mode wavefunction positioned mostly outside the con¬ 
fining potential, which reduces the impact of scattering. The mean free path rapidly grows 
with wavelength, in a direct analogy with guided optical waves in weakly guiding fiber de¬ 
signs, where weak confinement is employed to achieve exceptionally long mean free paths. 

Josephson junctions: Device overview 


We analyze five graphene Josephson junctions on hBN with widths ranging from W = 
800 — 1200 nm and lengths ranging from L = 250 — 350 nm (see Fig. Id for a labeled 
device schematic). Listed in Table SI are details on individual sample geometries. The 
small L/W aspect ratios place these devices are in the narrow junction limit, where the the 
critical current Jc can be approximated as a phase dependent summation over many parallel 
ID current channels (Equation 2 in the main text). Electrical measurements are conducted 
using standard Eockin techniques in a Eeiden Cryogenics Model Minikelvin 126-TOE dilu¬ 
tion refrigerator with a base temperature of 10 mK, well below the critical temperature of Al. 

Using a dry transfer method, graphene/hBN stacks are sequentially deposited on a 300 
nm thermally grown Si 02 layer, which covers a doped silicon substrate functioning as a 
global back gate. Graphene flakes are etched to the desired geometry using a 950 PMMA 
A4 polymer mask (~ 200 nm thick; spun at 4000 rpm) followed by an RIE 02 plasma etch. 
Titanium/aluminum (Ti/Al) superconducting electrodes are defined on selected flakes using 
electron beam (ebeam) lithography on a 950 PMMA A4 resist mask, followed by thermal 
evaporation and liftoff in acetone. Eor the titanium adhesion layer, we evaporate 10 nm at 
a rate of 0.3 Angstrom/s. This is followed by an evaporation of a 70 nm aluminum layer 
at a rate of 0.5 Angstrom/s at pressures in the low to mid 10“^ Torr range. Eor dual-gated 
bilayers, suspended top gates are fabricated using a standard PMMA/MMA/PMMA trilayer 
resist method which leaves a 200 nm air gap between the top gate and graphene. After us¬ 
ing ebeam lithography to define the gates, which employs position-dependent dosage, Cr/Au 
(3/425 nm) gates are deposited using thermal evaporation and liftoff in acetone. To remove 
processing residues and enhance quality, devices were current annealed in vacuum at dilution 
refrigerator temperatures. We note that edge currents were detected both in current-annealed 
and intrinsically high quality non-annealed devices; typically the appearance of edge currents 
coincided with the occurrence of Eabry-Perot interference in the ballistic transport regime. 
All five graphene Josephson junctions exhibit similar transport behavior. Additional data 
sets are provided in the Supplementary Eigures. 

Eourier method for extraction of supercurrent density distribution 


In a magnetic field B, the critical current Ic{B) through a Josephson junction equals the 
magnitude of the complex Eourier transform of the current density distribution J{x)'. 


UB) = |I,(B)| = 


J{x) exp( 27 r/(L -f lAi)Bx/^o)dx 


(28) 



where x is the dimension along the width of the superconducting contacts (labeled in Fig. 
Id), L is the distance between contacts, Iai is the magnetic penetration length (due to a 
finite London penetration depth in the superconductor and flux focusing), and $0 = hj^e 
is the flux quantum. Relevant in the narrow junction limit where current is only a function 
of one coordinate. Equation (28) provides a simple and concise description of our system. 
We employ Fourier techniques introduced by Dynes and Fulton to extract the real space 
current density distribution from the magnetic interference pattern By expressing the 

current density as J{x) = Js{x) + Ja{x), where Js{x) and Ja{x) are the symmetric and 
antisymmetric components, the complex critical current can be rewritten as: 

(29) 

/ CXD poo 

Js{x) cos{2t[{L + lAi)Bx/^o)dx + i / Ja{x)?,in{2T[{L + lAi)Bx/^Q)dx 

■00 J —00 

We calculate symmetric component of distribution, the relevant quantity for analyzing edge 
versus bulk behavior, as the antisymmetric component goes to zero in the middle of the sam¬ 
ple. For symmetric solutions, Ic{B) is purely real. To reconstruct Ic{B) from the measured 
critical current, the sign of Ic{B) is reversed for alternating lobes of the Fraunhofer inter¬ 
ference patterns. The extracted supercurrent distribution is expressed as an inverse Fourier 
transform: 


( 30 ) 


J s (x) 


Xc{B) exp(27ri(L + Iai)Bx/^ o)dB 


Because Ic{,B) is only nonzero over a rectangular window dictated by the finite scan range 
Bmin < B < Bmax, distribution extracted numerically is given by the convolution of J{x) 
with the sine function. To reduce artifacts due the convolution, we employ a raised cosine 
Alter to taper the window at the endpoints of the scan. Explicitly, 

( 31 ) Js{x) K, f Ic{B) cos^{7rB/2 Lb) exp{27ri{L + lAi)Bx/^o)dB 

^ ^min 

where n = 0.5 — 1 and Lb = {Bmax — Bmin)/2 is the magnetic field range of the scan. 


Gaussian fits to extract edge state widths 


To extract a length scale for the width of the edge currents near the Dirac point, we fit the 
experimental supercurrent density distribution Jc{x) to the Gaussian function 


(32) 




-f exp 


— {x + a)^ 
c 


where a determines the spatial peak offset, b determines peak height, and c determines 
peak width. For the data in Fig. IH, the fit parameters are a = 0.515, b = 8.8, and 
c = 0.017. The effective edge current width, given by the Gaussian full width at half maxi¬ 
mum xpwHM = 2.Vc ■ In 2, is 220 nm. 


Edge versus bulk amplitudes 
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To more quantitatively assess the evolution of edge and bulk currents with electronic car¬ 
rier density n, we plot line cuts of the individual contributions (see Fig. 2f and 3b). These 
are given by: 


(33) 




—x\Y+e\ 

E 

Xi — ^ VK 


Jc{,Xi,n) 


and 


Jr (n) = 


£2 


Xi = -£2 


Jc{.Xi,n) 

N2 


for a graphene flake whose full width spans from —xw to xw- the spatially aver¬ 

aged current amplitude over a small window of width ei from the edge of the flake. Similarly, 
jbui'^r) is the spatially averaged current amplitude over a strip of width 2^2 around the cen¬ 
ter of the flake. Ni = eijxstep and N 2 = £2 lx step, where Xgtep is the distance between data 
points (determined by the magnetic field range of the scan). For example, for the plots in 
Fig. 2F, Xw = 405 nm, ei = 29 nm, and £2 = 87 nm. 


Based on the edge versus bulk current profiles, one may infer whether edge doping is the 
dominant cause of edge currents in our devices. In the presence of edge doping, the edge ver¬ 
sus bulk contributions should be reversed for opposite polarities of bulk carriers (for example, 
edge dominated behavior at high densities on the electron side and bulk dominated behavior 
at high densities on the hole side), which is not consistent with the data. Bulk-dominated or 
flat distributions appear at both high electron and hole doping fairly consistently. As a sec¬ 
ond test, one can track the edge versus bulk contributions through the Dirac point to detect 
an offset in gate voltage between the charge neutrality point at the edge versus in the bulk. 
We did not detect positional density offset substantial enough to account for the large edge 
currents in these devices (Fig. 2F). 


Bayesian method for extraction of supercurrent density distribution 


The critical current as a function of the magnetic field, Ic{B), is related to the current 
density through the junction, Jc{x), as 


(34) 


Ic{B) = / dx Jc{x) exp {2TTixLBl^o) 


W 

2 


with L and W the length and width of the junction, and $0 = /i/2e the superconducting flux 
quantum. 


In the measured \Ic{B) \ all information about its complex phase is lost, making the prob¬ 
lem of determining the current density not have a unique solution. Using the method of 
Dynes and Fulton (DF), a unique solution can be found under the assumption of a symmetric 
current distribution, Jc{x) = Jc{—x). In practice however, disorder and inhomogeneities in 
the junction will lead to asymmetric current densities. Additionally, since experiments are 
performed over a finite range of magnetic fields, there is a cutoff in the current density reso¬ 
lution. Neither this finite resolution, nor experimental uncertainties are taken into account in 
the DF method, meaning it can only provide a qualitative estimate of Jc{x). 
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To gain a more quantitative understanding, we instead ask what is the distribution of Jc{x) 
whieh produees the same eritieal eurrent Ic{B). We answer this question by performing 
Bayesian inferenee to obtain the posterior distribution of the eurrent density, given the mea¬ 
sured eritieal eurrent. In our ease, Bayes’ rule reads: 


(35) 




P(|/e|; Jc)P(Jc) 


Here, P( Jc; l-^d) is the posterior distribution of the eurrent density, the quantity we want to 
ealeulate, while P{Jc) is its prior distribution. The likelihood funetion P(|/c|; Jc) indieates 
the eompatibility of the measured eritieal eurrent with a given eurrent density: 


(36) 


P(|/c|; Jc) = exp 


2^2 


where // is the eurrent obtained from Jc by using Eq. (34), Ic is the measured eurrent, and e 
is the measurement error. The faetor P(|/c|) is the same for all eurrent densities, meaning it 
does not enter into determining their relative probabilities. 


The experimental eurrent profiles are extraeted from seans of the differential resistanee as 
a funetion of DC eurrent bias and magnetie field, dV/B). Within the same sean, for 
some field values dV/dl has a elear maximum, while for others it monotonieally inereases 
towards its normal state value. We extraet the eritieal eurrent as the value Jdc at whieh the 
differential resistanee is a; x max dV/dl, ehoosing a value of a; < 1. This seleets points 
elose to the maxima at field values where they are well defined, and points elose to where the 
differential resistanee reaehes its normal state value otherwise. The uneertainty is obtained 
in the same fashion, by ehoosing a slightly smaller eutoff. 


We maximize the likelihood funetion using a Monte Carlo sampling algorithm [1]. To get 
a large resolution of the eurrent density without a signifieant inerease in the dimensionality 
of the sampling spaee, we expand Jc{x) as 


N 

(37) Jc{x) = An cos{27rnx/L) 

n=0 

and enforee Jc{x) > 0 for all x. The An coeffieients determine the shape of the distribution, 
whieh in Eq. (37) is assumed to be symmetrie, Jc{x) = Jc{—x). Using an asymmetrie form 
would typieally lead to a eritieal eurrent whieh shows node lifting - the minima of Ic{B) 
have nonzero values. While this feature is present in the measured eritieal eurrent, it ean be 
aeeounted for by faetors other than an asymmetrie eurrent distribution [2], sueh as relatively 
small aspeet ratios (~5), and a non-sinusoidal eurrent-phase relationship arising from a large 
junetion transpareney. Using a symmetrie Jc avoids this ambiguity, and has the additional 
advantage of providing a more direet eomparison between our method and that of Dynes and 
Eulton. 
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The likelihood function is maximized by allowing the coefficients to vary at each 
Monte Carlo step. As N is increased the posterior distribution of the current density widens, 
an indication of over-fitting. This increase in uncertainty serves as a criterion for choosing 
N, which for the typical dataset is between 4 and 8. The priors of An are set to the uniform 
distribution [— max(/c), max(/c)]. 

An example of our method is shown in Fig. S4, using N = 5. The current density is peaked 
at the edges of the sample, a feature also recovered in the DF approach. The corresponding 
critical current is in good agreement with the measured one, with the exception of the regions 
close to the nodes. Fig. S4 indicates that the supercurrent through the junction flows mainly 
along its edges. As a further test of the edge state contribution, we modify the functional form 
of the current density in Eq. (37), to explicitly allow for edge states. We add delta functions 
to the current density at the edges of the sample, Jdx) Jc{x) + dL6{x + W/2) + dji6{x — 
W/2), and estimate the contribution of edge states as the ratio of di + dji to the total current 
density As the carrier density approaches zero a significant fraction of the supercurrent 
is carried by the edge states, with (di + dn )/ ~ 0.45 (see Fig. S5). 
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Supplementary Figures 




Figure S 1. 


Fig. SI. 

The excess contribution to the spatially-resolved DOS near a line delta function poten¬ 
tial, AN{e,x) = N{e,x) — No{e) vs. distance from the delta function. Subtracted is the 
bulk contribution Nq given in Eq.(S5). The left panel shows the full excess contribution ob¬ 
tained from Eq.(S6), the right panel shows the contribution solely due to the guided modes, 
Eq.(S18). The two contributions are nearly identical, confirming that the peak in DOS can 
serve as a telltale of the guided modes.Parameter values used: A = —1.5hv, energies e = 
0.5£o 7 O.lso, where Sq = TTh^/nn^, uq = lO^^cm”^ (higher peaks correspond to higher 
energy e values). 
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Figure S2. 


Fig. S2. 

(A) Edge-dominated SQUID-like interferenee pattern at neutrality in deviee MLl (n = 
2.38 X 10® em“^; eolorseale is dV/dI{VL)). From Fig. 2A in main text. (B) Real-spaee 
image of eurrent flow eonfined to the boundaries, from data in part (A). (C) Conventional 
Fraunhofer pattern for uniform eurrent flow at high eleetron density (n = 7 x 10^^ em“^). 
From Fig. 2E in main text. (D) Real-spaee image of eurrent flow eonfined to the boundaries, 
from data in part (C). 
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Figure S3. 


Fig. S3. 

(A) Sequence of Fraunhofer measurements in bilayer device BL3 for the current maps in 
panels (B) and (C), shown in plots of dV/dliVL) as a function of magnetic field B (mT) and 
current bias Idc (nA). (B) Real space image of current flow J{x) as a function of carrier 
density on the hole side, showing edge currents near the Dirac point and a continuous evo¬ 
lution of bulk flow. (C) Individual line cuts of J{x) plotted from (B). This is the data set in 
Fig. 4A, plotted with a properly scaled vertical axis (supercurrent density, nAZ/im). 










15 



Figure S4. 


Fig. S4. 

Bayesian estimation of the supereurrent distribution. Posterior distribution of the eurrent 
density near the Dirae point in deviee BL3 (left panel), and eorresponding eritical eurrent 
(right panel). The values of Ic obtained from the posterior distribution (orange) are in good 
agreement with the measured eritieal eurrent (blue). 
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Figure S5. 


Fig. S5. 

Ratio of the supercurrent carried by the edge states as a function of carrier density in device 
BL3 over the density range of n ~ —1 to —2.9 x 10^^ cm“^. Each scan corresponds to a 
Fraunhofer pattern, with Fig. S4 showing the 8**^ scan. (Increasing scan number corresponds 
to decreasing carrier density.) 
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Supplementary Tables 


Device 

L (nm) 

W (nm) 

Aspect ratio, LAV 

Contact width (nm) 

BLl 

250 

1200 

0.208 

400 

MLl 

300 

1200 

0.25 

300 

BL2 

300 

800 

0.375 

400 

BL3 

350 

1200 

0.292 

600 

BL4 

250 

900 

0.278 

400 


Table SI. 

List of device dimensions for the graphene Josephson junctions studied in this work. L 
and W refer to junction length and width, respectively, as labeled in Fig. Id of the main 
text. Contact width refers to the size of the superconducting Ti/Al electrodes in the direc¬ 
tion perpendicular to W. BLx and MLx refer to bilayer and monolayer graphene devices, 
respectively. 
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